{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# -*- coding:utf-8 -*-\n",
    "# 数据预处理\n",
    "# 数据清洗：过滤掉非居民用电数据，过滤节假日的用电数据\n",
    "# 缺失值处理：补充缺失的数据\n",
    "# 三种方法：Lagrange插值法和Newton插值法以及Series自带的interpolate\n",
    "#1、Lagrange插值法和Newton插值法解决实际问题中关于只提供复杂的离散数据的函数求值问题，通过将所考察的函数简单化，构造关于离散数据实际函数f（x）的近似函数P（x），从而可以计算未知点出的函数值，是插值法的基本思路。\n",
    "#2、实际上Lagrange插值法和Newton插值法是同一种方法的两种变形，其构造拟合函数的思路是相同的，而实验中两个实际问题用两种算法计算出结果是相同的。\n",
    "#3、实验所得结果精确度并不高，一方面是因为所给数据较少，另一方面也是主要方面在Win32中C++中数据类型double精度只有7位，计算机在进行浮点运算时截断运算会导致误差。实际问题中，测量数据也可能导致误差。\n",
    "#4、在解决实际问题中，更多是利用精确且高效的计算机求解。所以解决问题时不仅要构造可求解的算法，更重要是构造合理的可以编写成程序由计算机求解的算法，而算法的优化不仅可以节省时间空间，更能得到更为精确有价值的结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 5),\n",
       " (0, 12),\n",
       " (1, 3),\n",
       " (1, 10),\n",
       " (1, 19),\n",
       " (1, 20),\n",
       " (2, 4),\n",
       " (2, 7),\n",
       " (2, 10),\n",
       " (2, 16)]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "inputfile = 'missing_data.xls' #输入数据\n",
    "data = pd.read_excel(inputfile, header=None)\n",
    "# 判断DataFrame中的空值的位置\n",
    "sij = []\n",
    "for j in range(len(data.columns)):\n",
    "    s1 = data[data.columns[j]].isnull()\n",
    "    for i in range(len(s1)):\n",
    "        if s1[i]==True:\n",
    "            sij.append((j,i))\n",
    "sij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>B</th>\n",
       "      <th>C</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>235.833300</td>\n",
       "      <td>324.034300</td>\n",
       "      <td>478.323100</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>236.270800</td>\n",
       "      <td>325.637900</td>\n",
       "      <td>515.456400</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>238.052100</td>\n",
       "      <td>328.089700</td>\n",
       "      <td>517.090900</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>235.906300</td>\n",
       "      <td>203.462116</td>\n",
       "      <td>514.890000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>236.760400</td>\n",
       "      <td>268.832400</td>\n",
       "      <td>493.352591</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>237.151181</td>\n",
       "      <td>404.048000</td>\n",
       "      <td>486.091200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>237.416700</td>\n",
       "      <td>391.265200</td>\n",
       "      <td>516.233000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>238.656300</td>\n",
       "      <td>380.824100</td>\n",
       "      <td>493.342382</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>237.604200</td>\n",
       "      <td>388.023000</td>\n",
       "      <td>435.350800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>238.031300</td>\n",
       "      <td>206.434900</td>\n",
       "      <td>487.675000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>235.072900</td>\n",
       "      <td>237.348072</td>\n",
       "      <td>609.193564</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>235.531300</td>\n",
       "      <td>400.078700</td>\n",
       "      <td>660.234700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>235.314951</td>\n",
       "      <td>411.206900</td>\n",
       "      <td>621.234600</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>234.468800</td>\n",
       "      <td>395.234300</td>\n",
       "      <td>611.340800</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>235.500000</td>\n",
       "      <td>344.822100</td>\n",
       "      <td>643.086300</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>235.635400</td>\n",
       "      <td>385.643200</td>\n",
       "      <td>642.348200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>234.552100</td>\n",
       "      <td>401.623400</td>\n",
       "      <td>618.197198</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>236.000000</td>\n",
       "      <td>409.648900</td>\n",
       "      <td>602.934700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>235.239600</td>\n",
       "      <td>416.879500</td>\n",
       "      <td>589.345700</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>235.489600</td>\n",
       "      <td>420.748600</td>\n",
       "      <td>556.345200</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>236.968800</td>\n",
       "      <td>408.963200</td>\n",
       "      <td>538.347000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "             A           B           C\n",
       "0   235.833300  324.034300  478.323100\n",
       "1   236.270800  325.637900  515.456400\n",
       "2   238.052100  328.089700  517.090900\n",
       "3   235.906300  203.462116  514.890000\n",
       "4   236.760400  268.832400  493.352591\n",
       "5   237.151181  404.048000  486.091200\n",
       "6   237.416700  391.265200  516.233000\n",
       "7   238.656300  380.824100  493.342382\n",
       "8   237.604200  388.023000  435.350800\n",
       "9   238.031300  206.434900  487.675000\n",
       "10  235.072900  237.348072  609.193564\n",
       "11  235.531300  400.078700  660.234700\n",
       "12  235.314951  411.206900  621.234600\n",
       "13  234.468800  395.234300  611.340800\n",
       "14  235.500000  344.822100  643.086300\n",
       "15  235.635400  385.643200  642.348200\n",
       "16  234.552100  401.623400  618.197198\n",
       "17  236.000000  409.648900  602.934700\n",
       "18  235.239600  416.879500  589.345700\n",
       "19  235.489600  420.748600  556.345200\n",
       "20  236.968800  408.963200  538.347000"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# <1> 拉格朗日插值法\n",
    "# 自定义列向量插值函数\n",
    "# s为列向量， n为被插值的位置， k为取前后的数据个数， 默认为5\n",
    "from scipy.interpolate import lagrange #导入拉格朗日插值函数\n",
    "data1 = pd.read_excel(inputfile, header=None,names=['A','B','C'])\n",
    "\n",
    "def plotinterplate_columns(s, n, k=5):\n",
    "    y = s[list(range(n-k,n) + list(range(n+1, n+1+k)))]\n",
    "    y = y[y.notnull()]#剔除空值\n",
    "    return lagrange(y.index, list(y))(n)#向位置n出插值并返回该插值结果\n",
    "\n",
    "# 逐个判断每列是否需要插值\n",
    "lagij = []\n",
    "for i in data1.columns:\n",
    "    for j in range(len(data1)):\n",
    "        if (data1[i].isnull())[j]:\n",
    "            data1[i][j] = plotinterplate_columns(data1[i],j)\n",
    "            lagij.append((i,j,data1[i][j]))\n",
    "            \n",
    "\n",
    "data1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data1.to_csv('lagrange.csv')\n",
    "from pandas import Series\n",
    "Series(lagij).to_csv('lagij.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# <2> 牛顿插值法\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\"\"\"\n",
    "@brief:   计算n阶差商 f[x0, x1, x2 ... xn] \n",
    "@param:   xi   所有插值节点的横坐标集合                                                        o\n",
    "@param:   fi   所有插值节点的纵坐标集合                                                      /   \\\n",
    "@return:  返回xi的i阶差商(i为xi长度减1)                                                     o     o\n",
    "@notice:  a. 必须确保xi与fi长度相等                                                        / \\   / \\\n",
    "          b. 由于用到了递归，所以留意不要爆栈了.                                           o   o o   o\n",
    "          c. 递归减递归(每层递归包含两个递归函数), 每层递归次数呈二次幂增长，总次数是一个满二叉树的所有节点数量(所以极易栈溢出)                                                                                     \n",
    "\"\"\"\n",
    "\n",
    "def get_order_diff_quot(xi = [], fi = []):\n",
    "    if len(xi) > 2 and len(fi) > 2:\n",
    "        return (get_order_diff_quot(xi[:len(xi) - 1], fi[:len(fi) - 1]) - get_order_diff_quot(xi[1:len(xi)], fi[1:len(fi)])) / float(xi[0] - xi[-1])\n",
    "    return (fi[0] - fi[1]) / float(xi[0] - xi[1])\n",
    "\"\"\"\n",
    "\n",
    "@brief:  获得Wi(x)函数;\n",
    "         Wi的含义举例 W1 = (x - x0); W2 = (x - x0)(x - x1); W3 = (x - x0)(x - x1)(x - x2)\n",
    "@param:  i  i阶(i次多项式)\n",
    "@param:  xi  所有插值节点的横坐标集合\n",
    "@return: 返回Wi(x)函数\n",
    "\"\"\"\n",
    "\n",
    "def get_Wi(i = 0, xi = []):\n",
    "    def Wi(x):\n",
    "        result = 1.0\n",
    "        for each in range(i):\n",
    "            result *= (x - xi[each])\n",
    "        return result\n",
    "    return Wi\n",
    "\"\"\"\n",
    "@brief: 获得牛顿插值函数\n",
    "\"\"\"\n",
    "\n",
    "def get_Newton_inter(xi = [], fi = []):\n",
    "    def Newton_inter(x):\n",
    "        result = fi[0]\n",
    "        for i in range(2, len(xi)):\n",
    "            result += (get_order_diff_quot(xi[:i], fi[:i]) * get_Wi(i-1, xi)(x))\n",
    "        return result\n",
    "    return Newton_inter\n",
    "# 自定义列向量插值函数\n",
    "# s为列向量， n为被插值的位置，k为取前后的数据个数， 默认为5\n",
    "def plotnewton_columns(s):\n",
    "    y = s[s.notnull()]#剔除空值\n",
    "    Nx= get_Newton_inter(y.index, list(y))         #向位置n出插值并返回该插值结果\n",
    "    return Nx\n",
    "\n",
    "import pandas as pd\n",
    "from scipy.optimize import newton #导入牛顿插值函数\n",
    "inputfile = 'missing_data.xls' #输入数据\n",
    "data2 = pd.read_excel(inputfile, header=None,names=['A','B','C'])\n",
    "\n",
    "newij = []\n",
    "for i in data2.columns:\n",
    "    Newton = plotnewton_columns(data2[i])\n",
    "    for j in range(len(data1)):\n",
    "        if (data2[i].isnull())[j]:\n",
    "            data2[i][j] = Newton(j)\n",
    "            newij.append((i,j,data2[i][j]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[('A', 5, 227.97718197331577),\n",
       " ('A', 12, 237.58023935829542),\n",
       " ('B', 3, 969.78391168067287),\n",
       " ('B', 10, 194.83076304210508),\n",
       " ('B', 19, 6553214.130850032),\n",
       " ('B', 20, 50150578.655558825),\n",
       " ('C', 4, 494.22872791561548),\n",
       " ('C', 7, 486.68478798271082),\n",
       " ('C', 10, 611.33674198392237),\n",
       " ('C', 16, 655.49360854653901)]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data2.to_csv('newton.csv')\n",
    "from pandas import Series\n",
    "Series(newij).to_csv('newij.csv')\n",
    "newij"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# <3> Series自带的iterplote\n",
    "data3=pd.read_excel('missing_data.xls',header=None,names=['A','B','C'])\n",
    "\n",
    "###利用Pandas中interpolate进行缺失值的补充\n",
    "data_out = data3.interpolate()\n",
    "\n",
    "df = data_out - data3.fillna(0)\n",
    "Serij = []\n",
    "for i in df.columns:\n",
    "    for j in range(len(df)):\n",
    "        if df[i][j] != 0:\n",
    "            Serij.append((i,j,df[i][j]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "data_out.to_csv('Series.csv')\n",
    "Series(Serij).to_csv('Serij.csv')\n",
    "N_L_S = pd.concat([data1,data2,data_out],axis=1)\n",
    "N_L_S.to_csv('nls.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## 比较三者的填充结果\n",
    "N = pd.DataFrame(newij,columns = ['col', 'row', 'value'])\n",
    "L = pd.DataFrame(lagij,columns = ['col', 'row', 'value'])\n",
    "S = pd.DataFrame(Serij,columns = ['col', 'row', 'value'])\n",
    "N.set_index(['col','row'],inplace = True)\n",
    "L.set_index(['col','row'],inplace = True)\n",
    "S.set_index(['col','row'],inplace = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "D:\\Anaconda2\\lib\\site-packages\\ipykernel_launcher.py:3: DeprecationWarning: \n",
      ".ix is deprecated. Please use\n",
      ".loc for label based indexing or\n",
      ".iloc for positional indexing\n",
      "\n",
      "See the documentation here:\n",
      "http://pandas.pydata.org/pandas-docs/stable/indexing.html#deprecate_ix\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style>\n",
       "    .dataframe thead tr:only-child th {\n",
       "        text-align: right;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: left;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th>Newton</th>\n",
       "      <th>lagrange</th>\n",
       "      <th>Series</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>col</th>\n",
       "      <th>row</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"2\" valign=\"top\">A</th>\n",
       "      <th>5</th>\n",
       "      <td>227.977</td>\n",
       "      <td>237.151</td>\n",
       "      <td>237.089</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>237.580</td>\n",
       "      <td>235.315</td>\n",
       "      <td>235.000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">B</th>\n",
       "      <th>3</th>\n",
       "      <td>969.784</td>\n",
       "      <td>203.462</td>\n",
       "      <td>298.461</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>194.831</td>\n",
       "      <td>237.348</td>\n",
       "      <td>303.257</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>6553214.131</td>\n",
       "      <td>420.749</td>\n",
       "      <td>416.880</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>50150578.656</td>\n",
       "      <td>408.963</td>\n",
       "      <td>416.880</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th rowspan=\"4\" valign=\"top\">C</th>\n",
       "      <th>4</th>\n",
       "      <td>494.229</td>\n",
       "      <td>493.353</td>\n",
       "      <td>500.491</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>486.685</td>\n",
       "      <td>493.342</td>\n",
       "      <td>475.792</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>611.337</td>\n",
       "      <td>609.194</td>\n",
       "      <td>573.955</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>655.494</td>\n",
       "      <td>618.197</td>\n",
       "      <td>622.641</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "               Newton lagrange   Series\n",
       "col row                                \n",
       "A   5         227.977  237.151  237.089\n",
       "    12        237.580  235.315  235.000\n",
       "B   3         969.784  203.462  298.461\n",
       "    10        194.831  237.348  303.257\n",
       "    19    6553214.131  420.749  416.880\n",
       "    20   50150578.656  408.963  416.880\n",
       "C   4         494.229  493.353  500.491\n",
       "    7         486.685  493.342  475.792\n",
       "    10        611.337  609.194  573.955\n",
       "    16        655.494  618.197  622.641"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compare = pd.merge(N,pd.merge(L,S,left_index = True, right_index = True),left_index = True, right_index = True)\n",
    "compare.rename(columns={'value':'Newton','value_x':'lagrange','value_y':'Series'},inplace= True)\n",
    "compare.ix[:,:] = compare.ix[:,:].applymap(lambda x: '%.3f' % x)  # 一旦执行这个操作，数据类型转成字符串型\n",
    "# compare['Newton'] = compare['Newton'].astype('float64') # 这一句不能要，否则又恢复原样了\n",
    "# compare = compare.round(3)# 此句话对科学计数法没起作用\n",
    "compare.to_csv('compare.csv')\n",
    "compare"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
